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We discuss an efficient approach to the calculation of the internal energy in numerical simulations 

of spin systems with long-range interactions. Although, since the introduction of the Luijten-Blote 
algorithm, Monte Carlo simulations of these systems no longer pose a fundamental problem, the 
energy calculation is still an 0(N ) problem for systems of size N. We show how this can be reduced 
to an 0(N log N) problem, with a break-even point that is already reached for very small systems. 
This allows the study of a variety of, until now hardly accessible, physical aspects of these systems. 
_^ ■ In particular, we combine the optimized energy calculation with histogram interpolation methods 

to investigate the specific heat of the Ising model and the first-order regime of the three-state Potts 
model with long-range interactions. 
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The numerical study of systems with long-range interactions is notoriously difficult, due to the large number of 
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interactions that has to be taken into account. Specifically, the number of operations to calculate the energy of a 
single particle scales as the total number of particles in the system, in contrast to the case of short-range interactions, 
where the corresponding number of operations is of order unity. This implies that Monte Carlo-based methods are 
generally restricted to very small system sizes, which are still hampered by strong finite-size effects. Some years 
ago, this problem was resolved for the case of O(n) spin systems with (ferromagnetic) long-range interactions, for 
which a dedicated cluster algorithm was developed jj| . Since the efficiency of this algorithm is independent of the 
rC number of interactions per spin, speed improvements of several orders of magnitude could be obtained compared to a 

r~^ ' conventional cluster algorithm. This speed-up pertains to the generation of independent configurations, for which the 
calculation of the energy is not required. Indeed, a variety of interesting physical results could be obtained by means 
of this method, see, e.g., Refs. ||-||]. Whenever one needs to sample the internal energy, however, the improvement is 
much less dramatic: the major remaining advantage is that one only has to calculate the energy for truly independent 
configurations, rather than in every Monte Carlo step. Whereas this still implies that one can study systems which 
are an order of magnitude larger than those that can be accessed via Metropolis- type simulations (cf. Ref . || ) , one 
is eventually limited by the fact that the total computing time scales quadratically with the system size. One major 
disadvantage of this inaccessibility of the energy is the fact that it is not possible to apply histogram interpolations 
in order to obtain information on thermodynamic quantities over a large parameter space ||. In this paper, we point 
out that, for systems with periodic boundary conditions, this problem can be circumvented by calculating the internal 
energy in momentum space. Thus, one can apply a Fast Fourier Transform (FFT), reducing the total computational 
effort to 0(N log N) for a system containing N spins. Indeed, this is a natural choice if one recognizes that the total 
energy is just given by a (discrete) convolution, which is one of the major applications of the FFT. Remarkably, the 
computational overhead entailed by the FFT turns out not to be a limiting factor: already for very small systems it 
is more efficient than a direct calculation of the energy. 

The remainder of this paper is organized as follows. First, we derive an expression for the energy in terms of the 
Fourier-transformed spin system and point out that several other observables can be obtained on the fly, at negligible 
additional cost. We also give a detailed comparison of our approach and the conventional method. Next, we illustrate 
our approach by means of several new physical results for one-dimensional systems with long-range interactions. We 
end with a summary of our results. 



II. ENERGY CALCULATION 



We will now first illustrate our approach for a (d = l)-dimcnsional system with an n-component order parameter, 
i.e., a generalized O(n) spin chain. This system is described by the Hamiltonian 
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W = -2^^ J(:E " 2/)S(x) ' S(y) ' (1) 
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where the spins S(x) are n-component unit vectors and N is the system size. The interaction J(x) is defined for all 
xeN. Under the condition that periodic boundary conditions are employed, the effective coupling J{x) between two 
spins is given by the sum over all periodic copies, 

oc 

j(x) = Yl J ( x + mN ) ( 2 ) 

m— — oo 

and hence has a period N. We set the self energy J(mN), which is just an additive constant in the total energy, equal 
to zero. Each component of the spin configuration S(a;) and the interaction J(x) can then be written as a Fourier 
sum 

1 N ^ 

fW = 7^ E h^ ikx/N , (3) 

* fc=0 

where the Fourier coefficients fk are obtained from the discrete Fourier transform of f(x). By means of the discrete 
convolution theorem J7J it is then straightforward to show that Eq. (|l|) can be written as 

H = -^fj^J k S k -S- k . (4) 

fe=0 

The essential step is now, that application of the Fast Fourier Transform J7||| reduces the computational effort for 
the calculation of the N Fourier coefficients from 0(N 2 ) to O(NlogN), thus, in principle, greatly speeding up the 
calculation of the total energy. The sum in (||) adds another O(N) operations, but this is compensated for by the fact 
that one typically also wants to calculate the magnetic susceptibility A^^ 1 1 ^2 x Zq S(a;)| 2 , which in the momentum- 
space representation is immediately given by |Sfc = o| 2 - For maximum efficiency, one has to restrict the system size to 
powers of 2. Naturally, the calculation of the coefficients J k has to be carried out only once. 

Even more can be gained, if one also desires to calculate the spin-spin correlation function g(r) — (S(0) • S(r)) = 
TV -1 J2 x =o ^( x ) ' &(. x + r )- The di scre te correlation theorem |7J states that the Fourier transform g k of g(r) is equal 
to A r_1 / 2 Sfc ■ S_fe, so that g k is obtained by N multiplications rather than another 0{N 2 ) operations in the real-space 
representation. 

All the above estimates are only measures for the complexity of the algorithm, which become valid for sufficiently 
large N. It remains to be seen whether the FFT-based approach is actually faster for the range of system sizes that 
can be accessed in present-day Monte Carlo simulations, which for lattice models go up to N ~ 10 5 — 10 6 . Figure p] 
compares the required CPU time per spin for the calculation of the internal energy via Eq. (|l|) and Eq. (^) , respectively, 
and the susceptibility. As expected, the former method scales asymptotically linearly with N. For the latter method, 
two estimates are given, which only differ in the choice of the FFT implementation. The slower results (op en squares) 
were obtained by means of the routines of Ref . H and the faster (triangles) by means of those of Ref . Il0[ . Although 
these two estimates differ by as much as a factor of 5, both of them outperform the conventional method already for 
N > 10; for N = 2 18 the improvement amounts to roughly four orders of magnitude. The initial downward trend in 
Fig. [j] is due to overhead being distributed over an increasing number of spins. Likewise can the "irregularities" in 
the FFT estimates be attributed to computational aspects. The slight deviations from linearity in the conventional 
estimates, however, are due to statistical inaccuracies in the timing: for N > O(10 4 ) this method becomes already 
prohibitively slow. Thus, we conclude that the method presented here provides a very efficient approach to energy 
calculations in lattice systems with long-range interactions; while there is still a weak system-size dependence in the 
required computational effort per spin, this no longer constitutes a bottleneck for practical applications. Note that 
higher-dimensional models can also be treated in this fashion. 



III. APPLICATIONS 
A. The Ising chain 

As a first example, we consider the Ising chain with algebraically decaying interactions, J(x) = J|x| _1_CT . The 
critical behavior for this system is essentially classical for < a < | and nonclassical for | < a < 1, see Ref. [ pj . 
Numerical results for the thermal exponent yt ]12| ] have indicated that the latter regime can be subdivided into two 
parts: y t > \ for \ < a < 0.65 and y t < \ for interactions that decay faster. This implies that the specific heat 
only diverges in a part of the nonclassical regime and should display a cusp- like singularity in the remaining part. 
By means of illustration, we have calculated the specific heat for a = 0.25 and a — 0.90. In both cases, we expect 
to find a function that does not diverge at the critical point, although the behavior should be qualitatively different. 
Simulations were carried out for N = 2 P , with 3 < p < 16, at a number of different couplings, for several times 10 6 
independent samples per system size. The full curves were determined by means of the multiple-histogram method [pi, 
where great care was taken to minimize systematic errors due to the histogram interpolation. 

Figure |] shows the specific heat C for a = 0.25, as a function of the reduced coupling K = J/(/sbT). It displays sev- 
eral close similarities to the specific heat of the mean-field model, including the build-up of a jump discontinuity at the 
critical point, the crossing of the finite-size curves in a single point at K c (up to corrections to scaling), and (not visible 
on this scale) an excess peak in the curves for finite systems, i.e., liniAr^oo C max (7V) ^ X\yo.kik c liniAr^oo C(K, N) [ p^[ . 
As shown in Ref. |2J], the location of the specific- heat maximum shifts as a function of system size, according to 

^max = K c + ai L~y" + a 2 L-^ + b^- 1 + • • • , (5) 

where y t * = ^ and the coefficients a*, 6$ are nonuniversal. A fit to this expression yielded yl = 0.51 (6) and K c = 
0.1147 (5), in good agreement with K c = 0.114142 (2) [|13|. The inset shows the peak height as a function of system 
size, strongly suggesting that the maximum is indeed finite in the thermodynamic limit. 

The case a = 0.90, shown in Fig. 0, clearly exhibits a distinctly different behavior. The specific heat is now nonzero 
in the thermodynamic limit, on either side of the critical point, and indeed displays the expected cusp-like singularity. 
The inset confirms that the maximum is convergent for N — > oo. Since y t is still sufficiently close to ^ i.e., the 
absolute value of the exponent a is sufficiently small, the location of the maximum cannot be distinguished from the 
critical point, unlike the case a = 1, where it is expected to occur at a coupling K < K c . 

B. The three-state Potts chain 

The ferromagnetic Potts model provides a particular generalization of the Ising model with respect to the number q 
of possible coexisting ordered phases. For q = 2 the Ising model is recovered; for q > 2 the ferromagnetic Potts model 
defines a genuine universality class distinct from the Ising and the more general O(n) universality class. The Potts 
model is of particular theoretical interest, because the phase transition it describes may be of first or second order 
depending on q and the spatial dimension d, even in the absence of symmetry-breaking fields. For nearest-neighbor 
interactions in d — 2 many properties of the Potts model are exactly known JL4||. In particular, if the model is in the 
first-order regime and q is sufficiently large, the asymptotic finite-size properties of the nearest-neighbor ferromagnetic 
Potts model have been established in a rigorous fashion |l5J. For long-range interactions, however, much less is known 
and currently available numerical data are limited to rather small systems p| . We demonstrate in the following, that 
also for the Potts model the cluster algorithm introduced in Ref. |1J can be combined with the FFT, allowing the 
numerical treatment of much larger systems. 

Again, we concentrate on the case of algebraically decaying interactions J(x) = |x| _1_ °'. The Hamiltonian of the 
ferromagnetic Potts chain with periodic boundary conditions can then be written in the same form as Eq. (|lj), where 
the Potts spins S(x) are unit vectors which mark the corners of a (hyper)tetrahedron in q — 1 dimensions. For the 
present case q = 3 we employ the complex notation 

S(x) -► S(x) € {l, e 2 "/ 3 , e 4 "/ 3 } , (6) 

i.e., S(x) ■ S(y) = di[S(x)S(y)*], where the asterisk denotes the complex conjugate. The spin representation of the 
Potts model given by Eq. (|l|) is equivalent to the standard Kronecker representation, but it has the advantage that the 
configurational energy is directly accessible by means of the FFT. According to mean-field theory the ferromagnetic 
Potts model should always show a first-order phase transition for q > 2. For our case of d = 1 and algebraically 
decaying interactions, one therefore expects the mean-field prediction to be correct for sufficiently small values a > 



of the decay exponent of the interaction, i.e., there should be a critical value <j c separating first- and second-order 
behavior. Mean-field theory provides an important guideline for the interpretation of our Monte Carlo data, so we 
briefly summarize the basic mean- field predictions. Following Ref. |Q we introduce the probability p K (x) that lattice 
site x is occupied by the Potts state re, 1 < re < q, and we define a homogeneous scalar order parameter s indicating 
a broken symmetry with respect to Potts state re = 1: 

i \ 1 + (l ~ l) s r \ 1 — s ,_> 

pi(x) = mi = , p K {X) = m K = , 2 < re < q . (7) 

For a given value of s the mean- field free-energy density /mf(s) (in units of UbT) is then obtained as 

/mf(s) = -K((l + a)s 2 + {[l + (q- l) s ] log[l + (q - l)s] + (q - 1)(1 - s) log(l - s)} /q , (8) 

where K — J / (k^T) denotes the reduced coupling and C( a ) is the Riemann zeta function. Note that the replacement 
K — ► (q — Y)q~ 1 K transforms Eq. m) from the spin representation into the Kronecker representation. The transition 
point K = K^-p from the disordered phase s = to the ordered (re = 1) phase s = s^j F follows from standard 
mean- field arguments |f4[ : 

K t MF ai + a)= { -^-^-\og( q -l) : 4 F = i^. (9) 

q(q-2) q-1 

According to Eqs. ^h and (M) the distribution function P(mi) for a finite system displays three maxima near the 
transition temperature T/^p = J / (k^K^): one at m\ = 1/q for the disordered phase, one at mi = (q — \)/q for the 
ordered (n = 1) phase, and one at mi = l/[q(q— 1)] for the ordered phases with respect to the remaining Potts states 
(re > 2). Note that all ordered phases appear with equal probability in the course of the simulation. For a < 0.4 and 
in our case of q — 3 these three peaks in P(m\) are indeed located very close to their mean- field positions. For a = 0.6 
the peaks are still clearly separated, but they occur at positions shifted with respect to the mean-field predictions and 
for a > 0.7 the peaks start to overlap strongly and can only be identified for very large systems (see below). 

Although the algorithm introduced in Ref. ]lj is by far the most efficient one for the simulation with spin systems 
with long-range interactions, it is not able to deal with first-order phase transitions beyond a certain system size. The 
reason is that like the Metropolis algorithm the Wolff cluster algorithm encounters an activation barrier between states 
with and without long-range order, which is set by the energy-density gap between the disordered and the ordered 
phase. For a given size of the gap the tunneling time between disordered and ordered phases, and therefore the 
required sampling time, increases exponentially with the system size so that the attainable system size N is severely 
limited. This tunneling problem can be solved by employing the well-established ideas of multicanonical sampling p6[ , 
however, the generalization of the cluster algorithm O to an efficient multicanonical algorithm is beyond the scope 
of the present article. 

The data we present in the following have been obtained from histograms of the energy taken at several temperatures. 
The data are again conveniently analyzed by the optimized multiple- histogram method B. For the values a = 0.2 
and a = 0.4 the Potts chain undergoes a strong first-order phase transition p| which limits the chain length to 
TV = 2 13 in the former case and to N = 2 14 in the latter case. We reinvestigate chains from N = 2 10 spins to the 
respective maximum chain length by taking a few times f 6 independent samples for each system size and temperature, 



where a comparison with the finite-size theory of Borgs et al. |15| turns out to be very instructive. We restrict the 
detailed presentation to the case a — 0.4; the data for a = 0.2 are qualitatively very similar. Near the transition 
temperature the energy distribution function P(E) displays two peaks characterizing the ordered and the disordered 
phase, respectively |L7J. In Ref. fa] the temperature of equal peak height is taken as an estimate for the transition 
temperature on the finite system, where the leading finite-size corrections are of the order 0(1/N). If the systems 
are large enough, i.e., the peaks in the energy distribution are well separated, the ratio W /Wd of the weight of the 
ordered phase W and the weight of the disordered phase Wd provides a far more convenient indicator of the transition 
temperature, because the associated finite-size corrections decay exponentially with TV |lq , |l7| ]. Our result for long- 
range interactions is shown in Fig. ^. The transition temperature is marked by the intersection of W /Wd as a function 
of temperature for the three largest systems. For TV = 2 11 the peaks in the energy distribution are not well separated 
so that W /Wd is not well defined in this case. For TV > 2 12 the curves meet at W /Wd = 1-67 (2) as shown by the 
solid line. For a — 0.2 we find a corresponding intersection at W /Wd = 1-25 (2). For nearest-neighbor interactions 
in d > 2 and sufficiently large q the value W /Wd = q is expected to indicate the transition temperature |f^ , |i7| . 
Surprisingly, we find a much smaller value here which appears to increase with a . From Ref. J15| one furthermore 
expects that the curves displaying the energy density for different system sizes as a function of temperature exhibit 
an intersection close to the transition temperature, where the deviations are predicted to be exponentially small in TV. 
In Fig. || this situation is shown for long-range interactions with a = 0.4. A corresponding result has been obtained 



for a — 0.2. The energy densities intersect near the transition temperature found in Fig. [|, where the shifts between 
mutual intersections seem to be compatible with exponentially small finite-size effects. Still too few data are available 
for a quantitative analysis of these shifts, but finite-size effects of the order 1/N can be ruled out. The fourth-order 
energy cumulant t/4 defined by 

U 4 = (H 4 )/(H 2 ) 2 (10) 

is shown in Fig. O for different system sizes as a function of temperature, where Ti is the Hamiltonian given by Eq. (Fy). 
These cumulants should also intersect at the transition temperature in the limit TV — ► 00. The data displayed in 
Fig. ra show the expected tendency, but the finite-size corrections are much larger than those for the weight ratio or 
the energy density (see Figs. [|and g, respectively). For a = 0.2a corresponding result has been found. The systematic 
shift of the intersections of U4, for different system sizes is compatible with a 1/TV behavior as anticipated in Ref. Eq ] 
for nearest-neighbor interactions, but the present amount of data is too limited to give reliable quantitative evidence 
for this behavior. 

For nearest-neighbor interactions and periodic boundary conditions the energy density asymptotically obeys the 
scaling law (cf. Eq. (1) of Ref. @) 



ip/a t\ E d + E Ed — E 
E(0, L) ~ tanh 



*^v-w+^ 



(11) 



where (3 = 1/ (k^T) is the inverse temperature and f3t is the transition point. It is instructive to compare the scaling 
form given by Eq. ( |ll| ) with the data displayed in Fig. [|. The energies of the disordered phase Ed and the ordered 
phase E can be read off from the positions of the two maxima of the energy distribution function. It turns out 
that the data in Fig. g and their counterpart for a = 0.2 are in fact consistent with Eq. (|ll]) within the error 
bars, provided the number of states q on the right-hand side of Eq. (0) is replaced by the effective value q e s(cr — 
0.4) = W /Wd\p=p t = 1-67 measured in Fig. [| at the transition temperature or its counterpart q c ff(a = 0.2) = 1.25, 
respectively. The quantitative comparison of our data for a = 0.2, TV = 2 13 and a = 0.4, N = 2 14 with Eq. (O) is 
shown in Fig. 0. Within the statistical errors, the agreement is excellent except for larger values of the scaling variable 
(/3 — j3 t )N . These deviations are due to the fact that Eq. ( |PT| ) only holds asymptotically for sufficiently large systems. 
For finite systems additional finite-size corrections enter through the residual TV-dependence of Eg, E , and q c g(o~) 
which appear as parameters in Eq. (|ll|). Figure f 7 ] demonstrates that the finite-size effects in the three-state Potts 
chain with periodic boundary conditions and long-range interactions can be interpreted in terms of the Borgs-Kotecky 
theory |l|] for the nearest-neighbor Potts model in higher dimensions for an effective number of states <7eff(c). The 
physical meaning of q c g (a), however, remains unclear. The proof of Eq. (0) also requires the assumption that q is 
sufficiently large |15| , so q = 3 may not be sufficient for a quantitative comparison. On the other hand, numerical 
investigations have shown that the q > 5 nearest-neighbor Potts model in d = 2 fllTJ and the q — 3 nearest-neighbor 
Potts model in d = 3 [lj] follow the theory of Ref. jl5| very closely despite the small values of q. Further analytical 
and numerical studies are required to settle this question. 

We close our discussion of the three-state Potts chain with a brief summary of our results for a = 0.6, 0.7, and 0.75, 
which have been studied with reduced statistics (10 5 independent samples for each system size and temperature). The 
value a = 0.6 is still located in the first-order regime g, but in order to obtain a well-defined weight ratio W /Wd 
system sizes of N > 2 16 Potts spins are required, although the maxima in the energy distribution are well separated 
already for N > 2 14 . We have performed simulations for N = 2 14 up to TV = 2 17 at four to six temperatures for each 
system size. Even for TV = 2 17 the finite-size corrections are too large to identify an intersection of the energy densities 
as accurately as displayed in Fig. [|. Data acquisition for TV > 2 17 is strongly hampered by the energy gap so that we 
refrain from discussing our data for a — 0.6 in any more detail. In contrast to Ref. H the value a = 0.7 of the decay 
exponent can undoubtedly be identified as a member of the first-order regime. For system sizes TV > 2 16 the energy 
distribution function displays the typical double-peak structure which becomes sharper as the system size is increased 



at fixed temperature. We illustrate this in Fig. 
close to the transition point. The same analysis 



, where the data for P(E) are shown at T/T^p = 0.8095, which is 
ras been repeated for a = 0.75 and system sizes up to TV = 2 19 spins. 
Although P(E) also develops a plateau similar to the one displayed in Fig. for TV = 2 15 , no double-peak structure 
could be resolved up to TV = 2 19 so that a — 0.75 may already belong to the second-order regime of the three-state 
Potts chain with long-range interactions. However, the detection of a double-peak structure in P{E) for a given value 
of a is essentially a matter of attainable system size (see Fig. ||), so a c > 0.7 is the only safe conclusion here. 

IV. CONCLUSIONS 

The combination of the recently developed cluster algorithm M for systems with long-range interactions with the 
Fast Fourier Transform for the calculation of the configurational energy leads to a Monte-Carlo algorithm with a 



very high efficiency. In particular, the FFT allows to extend the attainable system sizes by two orders of magnitude 
in comparison with other approaches (cf. Ref. H). Histogram interpolation methods then allow the investigation of 
thermodynamic properties of these systems with unprecedented resolution. By construction the algorithm can only 
deal with first-order phase transitions up to a limited system size. In order to avoid this limitation the algorithm must 
be generalized to include multicanonical sampling. Here, we have demonstrated the potential of the algorithm for the 
Ising chain and the three-state Potts chain with algebraically decaying interactions. For completeness, it is mentioned 
that also for system sizes that are not integer powers of two a considerable gain can be obtained by performing the 
discrete Fourier Transform via, e.g., a prime- factor algorithm. 

For the Ising chain we have investigated the finite-size behavior of the specific heat in the classical regime for 
a = 0.25 and in the nonclassical regime for a — 0.9. In the former case the specific heat behaves essentially mean-field 
like, i.e., the expected discontinuity in the specific heat at the critical temperature in the thermodynamic limit builds 
up as the system size is increased. On the other hand, the choice a = 0.9 is expected to yield a negative specific-heat 
exponent, i.e., a cusp singularity should appear with increasing system size. Our numerical data confirm this behavior 
as well and clearly show the different shapes of the specific-heat curves in the two cases. 

The three-state Potts chain is expected to show a first-order phase transition for a < a c , where our results indicate 
that <7 C > 0.7. For a — 0.2 and a = 0.4, for which the q — 3 Potts chain displays a strong first-order phase 
transition, our data confirm the Borgs-Kotecky scenario of the first-order transition in Potts models with nearest- 
neighbor interactions in higher dimensions, provided the number q of states is replaced by the effective number of 
states q c ff(a) = W /Wd\p=p t < q which also enters the finite-size scaling form of the energy density near the transition 
temperature. For a = 0.6 the same behavior can be confirmed only on a semi-quantitative level, because much larger 
systems must be investigated in order to obtain sufficient resolution. The mechanism that leads to the reduction of 
the effective number of states and the physical interpretation of q c g (a) are not known. 
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log io ^ 
FIG. 1. Required CPU time (in microseconds on a Pentium-II 400 MHz) per spin for the calculation of the energy and 
the susceptibility, as a function of system size, for both the conventional method and two FFT-based methods. For a further 
discussion see the text. 
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FIG. 2. Specific heat of the Ising chain with interactions decaying as r~ ' , for system sizes 8 < N < 65536. One observes 
the appearance of a mean-field-like discontinuity at the critical coupling. The inset shows the peak height as a function of 
system size. 
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FIG. 3. Specific heat of the Ising chain with interactions decaying as r~ 1,90 , for system sizes 8 < N < 65536. With 
increasing N, a cusp-like singularity appears. The vertical line indicates the critical coupling. The inset shows the peak height 
as a function of system size. 
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FIG. 4. Ratio W /Wd of the weights of the ordered and the disordered phase as a function of temperature for the system 
sizes N = 2 11 (o), 2 12 (+), 2 13 (□), and 2 14 (x). The temperature is measured in units of the mean-field transition temperature 
Tmf- Statistical errors (one standard deviation, not shown) are smaller than the symbol sizes. The horizontal line marks the 
value of Wo/Wd at the intersection of the three largest systems (see main text). 
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FIG. 5. Energy E/(JN) per spin in units of the coupling constant J as a function of temperature for the system sizes 
N — 2 11 (o), 2 12 (+), 2 13 (□), and 2 14 (x). The temperature is measured in units of the mean-field transition temperature 
TjJif- Statistical errors (one standard deviation) are only shown when they exceed the symbol sizes. Within the statistical error 
the curves for the three largest systems intersect at the same temperature as W /Wd in Fig. m 
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FIG. 6. Fourth-order energy cumulant XJa as a function of temperature for the system sizes N = 2 1 (o), 2 1 (+), 2 (□), 
and 2 14 (x). The temperature is measured in units of the mean-field transition temperature Tm F . Statistical errors (one 
standard deviation) are only shown when they exceed the symbol sizes. The curves do not have a common intersection within 
the displayed temperature range indicating much larger finite-size effects than in Figs. and n. 
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FIG. 7. Energy density for a = 0.2, N = 2 13 (o) and for a 



2 

0.4, TV = 2 1 



(+) as functions of the scaling variable (/3 t — /3)iV 



in comparison with Eq. ( |ll| ) for g — » g c ff(cr = 0.2) = 1.25 (solid line) and q — » <j e ff(c = 0.4) = 1.67 (dashed line). Error bars 
on the numerical data (one standard deviation) are only shown when they exceed the symbol sizes. The inverse temperature 
/3 = l/(k B T) is given in units of 1/(/»b3mf). 
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FIG. 8. Energy distribution function P(E) for a = 0.7, T/T^ F = 0.8095, and system sizes N = 2 15 (o) and AT = 2 17 (- 
The typical double-peak structure remains invisible for TV < 2 16 spins but sharpens if N is increased at fixed T. 
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